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ABSTRACT 

We present the analysis of a suite of simulations run with different particle-and grid- 
based cosmological hydrodynamical codes and compare them with observational data 
of the Milky Way. This is the first study to make comparisons of properties of galaxies 
simulated with particle and grid-based codes. Our analysis indicates that there is 
broad agreement between these different modelling techniques. We study the velocity 
dispersion - age relation for disc stars at z = and find that fo ur of the simulations 
are more consistent with observations bv iHolmberg et all (|2008[ ) in which the stellar 
disc appears to und ergo continual/secu l ar hea ting. Two other simulations are in better 
agreement with the lQuillen fc Garnettl ()2001[ ) observations that suggest a "saturation" 
in the heating profile for young stars in the disc. None of the simulations have thin 
discs as old as that of the Milky Way. We also analyse the kinematics of disc stars at 
the time of their birth for different epochs in the galaxies' evolution and find that in 
some simulations old stars are born cold within the disc and are subsequently heated, 
while other simulations possess old stellar populations which are born relatively hot. 
The models which are in better agreement with observations of the Milky Way's stellar 
disc undergo significantly lower minor- merger/assembly activity after the last major 
merger - i.e. once the disc has formed. All of the simulations are significantly "hotter" 
than the Milky Way disc; on top of the effects of mergers, we find a "floor" in the 
dispersion that is related to the underlying treatment of the heating and cooling of 
the interstellar medium, and the low density threshold which such codes use for star 
formation. This finding has important implications for all studies of disc heating that 
use hydrodynamical codes. 

Key words: galaxies: formation — galaxies: evolution — Galaxy: thick and thin disc — 
methods: N-body simulations 



1 INTRODUCTION 

One of the major outstanding "grand challenges" facing as- 
trophysics for the coming decade is the unravelling of the 
underlying physics governing the formation and evolution 
of disc galaxies such as our own Milky Way. A principal 
difficulty resides in trying to accommodate the early col- 
lapse and violent merging history intrinsic to the canonical 
framework of "hierarchical assembly" of galactic structure 



with the apparent stability of what should be fairly fragile 
thin galactic discs. 

High Performance Computing (HPC) simulations 
of gravitational A^-body and hydrodynamical physics 
have become a primary tool with which to m odel galaxy 
forma tion 

Jl992h: ISummers et aL , 

ISteinmetz fc Mueller! (|l994l V ISommer-Larsen et all (2003 



in a cosmologica l context (e.g Katz et all 



lll993h: iNavarro fc White! | l99 
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lAbadi et all (12003d); IRobertson et all (12004); iBailin et all d2004h and the follow u p stu d y by IHolmberg et all 



Okamoto et 




(120051); iGovernato et all (j2007j) 
ISanchez-Blazauez et al 

simulations model the for- 
mation and evolution of disc galaxies within a Universe 
dominated by a Cold Dark Matter (CDM) component and 
a cosmological constant (A) . While powerful, the techniques 
employed are not without their problems; for example, 
the loss of angular momentum in the luminous component 
of disc galaxies is one of the major problems in most 
of the aforementioned cosmological simulations. In these 
simulations, gas cools efficiently via radiative processes, 
causing baryons to collapse rapidly during the earliest 
phases of the hierarchical clustering process. The luminous 
component ends up transferring angular momentum to the 
dark matter halo making the luminous component deficient 
in angular momentum. This is often referred to as th e 
"angular momentum probl em" (Nav arro fe Benj l|l99lf ): 
ISteinmetz fc Navarro! ((2002)). As a result, these simula- 
tions typically produce galaxies with an ove rly-dominant 
spheroi d component and an overly small disc (jAbadi et al.l 
l|2003al ) ; IScannapieco et al] (120091)') . in d i sagree ment with 
observations of disc galaxies ( Brook et al.l (|2004h ). 



2007), ISoubiran fc Girardl |200 



Soubiran et al 



2009); iQuillen fc Garnettl (|200lF 7 and iDehnen fc Binnev 



Another challenge facing disc galaxy formation in the 
ACDM scenario is the old age of the Milky Way's thin 
disc. This seems at odds with the heating that one ex- 
pects from merging and accretion events within a ACDM 
paradigm. Indeed, several studies of isolated discs being 
bombarded by satellites have shown that one would expect 
that the disc would b e destroyed, or at le a st severely heated, 
by acc r etion events dQuinn et al l (1 19931); iKazantzidis et al.l 
l|2008t ); IKazantzidis et all <j2009t ): iRead et all (|2008l )). Two 
rece nt studies have incl uded gas in the main disc, with 
one dMoster et all (|2010D ) finding a significant decrease in 
heating, by 25% - 40% for g as fractions o f 20% and 40% 
respectively, with the other (|Purcell et al.l (|2010l )) finding 
that the effects of gas are somewhat less dramatic. What 
is clear is that all studies which use contrived initial con- 
ditions which are bombarded with satellites are necessar- 
ily restricted in their application, both for the disc and 
satellites. For example, how best to assign an appropri- 
ate velocity dispersion, mass, and scalelength of the Milky 
Way disc at redshift of two, say? Were the stars already 
kinematically "hot" in this early disc, or have they been 
heated subsequently? Idealised studies with pre-formed discs 
can be powerful, but they do not address directly the is- 
sues pertaining to disc formation and how this relates to 
merger events as they occur within a hierarchical cosmol- 
ogy. Suffice to say that the exist ence of thin di s cs rem ains a 
challenge for ACDM cosmology. IStewart et al.l (|2009l ) argue 
that gas rich-mergers can explain the number of low mass 
galaxies on the blue sequence and mass-morphology rela- 
tion, but their analysis is not able to address the issue of 
the thinness of t he discs which surv i ve mergers. In fact, i t 
has been shown dBrook et alJ j2004h;ISprmgel et al l (12005 " 



IRobertson fc Kravtsovf ( 2008); IGovernato fc Broo: 



that gas rich mergers in simulations result in hot thick discs, 
with thin discs forming in the subsequent quiescent period. 

Observations of the kinematics of disc stars of our 
Galaxy have been carried out throughout the years in or- 
der to understand the mechanisms gove rning the forma- 
tion of the disc. These studies include INordstrom et all 



2008 



1998 



These studies however, have provided different pictures 
of the relationship between the ages of disc stars 
and the ir velocity disper si ons ( the age-dispersion re- 



lation) . 



Quillen fc Garnettl j|200ll ), using the data of 
lEdvardsson et aL I (|l993l ). found that vertical disc heating for 
the Milky Way saturates at a w ~ 20 kms -1 , with the value 
of dispersion virtually constant for stars of ages between ~2 
and ~9 Gyr. A discrete jump is apparent for stars with ages 
>9 Gyr which is generally interpreted to be the signature 
of the thick disc. Thus, this study supports the notion of a 
thick disc as a separate component to the thin disc, and sug- 
gests dif ferent formation scenari os for each component. By 
contrast , INords trom et all l|2004l ) and the follow up study of 
IHolmberg et al.l (|2007l ) advocate a picture in which the disc 
has undergone continual heating over the past ~10 Gyr. It 
is not clear from these later studies whether a thick disc 
should be considered as a separate component: firstly, the 
selection is biased toward thin disc stars, and secondly, if 
the thick disc is a separate component, it is possible that 
their result is driven by increasing contamination of their 
samp le by thick disc stars as older and older stars are exam- 
ined l|Navarro et al.l (|2010l )). One of the main points of con- 
tention in such studies, and a possible explanation for the 
different findings, is the difficu lt y in determining t he ag es 
of stars ( Anguiano et all (12009 ): lAumer fc Binnevl (|2009l) ). 
Further, Seabroke"'fc***Glmiorel ( 2007 ) showed that a power 
law fit as suggested in the Geneva-Copenhagen studies is 
statistically similar to disc heating models which saturate 
after ~4.5 Gyr, and is consistent with a minimal increase of 
a m for old stars. We will compare our simulations to both 
these data sets. 

The degree of heating of thin disc stars is cer- 
tainly complicated by an old, hot "thick disc" of stars. 
Since the disco very of the thick disc component of the 
Milky Way by iGilmore fc R"eldl j| 19831 ). several analyses 
of ages, abundances, and kinematics have indicated that 
the thick and thin discs are two disti n ct components 
e.glReid fc Maiewskil (ll993l) ; lNissenl (1 19951) ;IChiba fc Beersi 



20001); iBensbv et al l fcoOTl ): iBeers et al] (|2009l ) although 



see llvezicTet'liL l|2008l ) for an opposing view), where a thin 
disc whose stars have formed continuously over ~ 9 Gyr 
is superimposed on an old thick disc. The thick disc has 
also been shown to be a common co mponent in most, and 
possib ly all, observed spiral galaxies (|Yoachim fc Dalcantonl 
(2008)), at least in the sense that the light distributions 
are better fit by two functions rather than one. Recent ob- 
servations carried out measuring the stellar population of 
the Milky Way's thick disc give sc a leheig h ts that range 
betwe en 500 - 1100 pc l|juric et all l|200Sl) : ICarollo et all 
(2010)) compared to the thin disc that h as a measured scale - 
height that ranges b etween 200 - 400 pc l|juric et al.l (2008); 
ICarollo et al] ([2010)). The rotational lag o f the thick disc 
of the Milky Way ra nges from 20 kms" 1 l|Chiba fc Beersi 
l|2000l) ) to 50 kms -1 l|Soubiran et all l|2008l) ). The velocity 
ellipsoid (i.e., the velocity dispersions in the u,v,w refer- 
ence frame) of the thick disc has been quoted as being from 
(cr u , cr v , a w ) = ( 46±4 , 50 ± 4, 35 ± 3) kms~ 1 , as measured by 
IChiba fc Beersi (|2000l ). to (a u ,a v ,a w ) = (63± 6, 39±4, 39±4) 
kms -1 , as measured bv lSoubiran et ah! (2008). The thin disc 
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dominates the thick disc in the local region by a factor of 10:1 
in terms of stellar mass, but difficulty in determining scale- 
lengths of the two means that comparing their total masses 
is highly unce rtain, with e stima tes ranging from mass ratios 
of 10:1 to 3:l |Juric et ail (|2008l ). 

For this study we focus on the kinematics and heating 
of all stars within the disc region of our simulations, 
without any a priori distinction between the thin and 
thick discs. We will explore and discuss occasions where 
two components arise, in the hope of shedding light 
on the various scenarios which have been postulated to 
explain the formation of the thick and thin discs. One 
model suggests that a previously existing thin disc was 
kinematically heated to give rise to the thick disc. This 
vertical heating of the thin d isc st ars could be rapid, 
due to mergers (|Quinn et al.l |l993)), or more gradual, 
as a result of secular processes such as giant molecular 
clouds, spiral arms or the presence of a b ar prov i ding t he 
necessary kinematic "kick" to the stars l|Larsonl l| 19761) ). 
An alternative model suggests that the thick disc formed 
during the violent relaxation of the galactic potential prior 
to the formation of the thin disc, where star formation was 
triggered by accreti on of gas du r ing m ajor merger events 
at an early epoch ferook et al.l ((200J)) . This raises the 
question of whether the age-dispersion relation is at least 
in part due to the earlier discs being, in general, hotter 
than the later discs - i.e. the old stars were born hotter 
than the younger stars, possibly related to earlier mergers 
in ACDM having higher mass ratios between the mass of 
the central gal axy to the ac c reted satellite in general than 
later mergers ((Brook et alJ (120051 )). Recent observations 
suggesting that high-re ds hift d isc s are relatively thick 
Dalcant on fc Bernstein! (120021); van Starkenburg et al.l 



20081 ): iLemoine-Busserolle fc Lamareillel (|20icl ))" jossibly 
provide support for this scenario. A further alternative 
model suggest that the thick disc stars are actually accreted 
from sat ellite galaxie s durin g the hierarchical assembly 
process (|Abadi et all (|2003al )). Recently, models of discs 
which are entirely isolated from the satellite bombardment 
that is predicted in ACDM have shown that migration of 
stars can naturally lead to combinations of age, metallicity 
and dispersions which are consis tent with observations of 
thin and thick disc populatio ns (lLoebman et al. I 20ld) as 
well as an analytic model by ISchonrich fc Binnev (2009)). 
These models beg the question of whether heating is 
required at all, yet to be fair, they have not been in tegrated 
within a fully cosmological paradigm. Finally, iKroupal 
(2002) suggested that massive star clusters formed at high 
redshift dissolve at later times to form the thick disc, a 
theor y given support by the clumpy nature of high reds hift 
discs (jElmegreenl (|2007l ): lElmegreen fc Elmegreeni l|2006h ). 

We aim to provide further insight into the causes of the 
age-velocity dispersion relation of all disc stars (thick and 
thin) by using a suite of HPC simulations of disc galaxies, 
each run with a different hydrodynamical code, different ini- 
tial conditions, different resolution, and different assembly 
histories, to sample a wide range of disc galaxy formation 
pathways. In particular, we determine whether disc stars are 
born hot or cold in early times, and the degree in which they 
subsequently heat vertically. Attention is given to any con- 
nection between heating rates and accretion histories. The 
issue of the effects of numerical heating is important in all 



studies of disc heating. We show that numerical heating is 
not causing the measured disc heating. We highlight the 
role of the implementation of star formation recipes, and 
the modelling of the interstellar medium in which stars are 
formed, as determinants of a dispersion "floor" for the sim- 
ulations. 

In § 2, we briefly describe the main characteristics of 
each of the different codes used to produce our compilation 
of simulations. In § 3, we present our main results which 
include the dispersion of all stars as a function of time for the 
final galaxies, comparing the simulations with observations 
of the Milky Way, studying the kinematics of young stars 
at different epochs, as well as following the heating of stars 
born at early epochs; we then relate the heating with merger 
processes within the simulations. In § 4, we examine isolated 
simulations, and show that heating does not occur in the 
absence of a cosmological environment, ruling out numerical 
effects as the primary agent driving our results. We present 
our conclusions in 5 5. 



2 THE SIMULATIONS 

We analyse seven cosmological disc simulations run with dif- 
ferent iV-body hydrodynamical galaxy formation codes. In 
this section we provide a summary of the main details of 
each code. For full details, references are provided. 

Two of the si mulation s we a nalysed in this study are 
run with RAMSES l|Tevssieri (|2002l )). which models the gas 
hydrodynamics using an adaptive mesh refinement (AMR) 
scheme, while the other codes use a smoothed particle 
hydrodynamics (SPH) approach^ Using examples drawn 
from these different fundamental approaches should provide 
greater confidence in the robustness of our results. To the 
best of our knowledge, our study is the first to compare 
properties of simulated disc galaxies formed using these two 
commonly-adopted methodologies. 

All five simulation^] have cosmological initial condi- 
tions where small scale structures merge to form increasingly 
larger objects in the Universe as part of the so-called hier- 
archical framework. They all have a similar Milky Way-type 
mass halo and all the different codes self-consistently include 
the primary physical processes needed to model galaxy for- 
mation and evolution. These consist of the effects of grav- 
ity, hydrodynamic pressure and shocks, star formation and 
feedback, radiative cooling, and a photoionising UV back- 
ground. They all adopt a type of Schmidt law for converting 
gas particles into stars, where the star formation rate is pro- 
portional to the gas density to some power. 

The main difference between the simulations is that 
they have different initial conditions, and hence merger and 
assembly histories. They also adopt different recipes for feed- 
back from Type II supernovae (SNell). Various methods 
have been suggested to incorporate the supernovae feedback 



1 iKnebd J2005f) provides an excellent primer to the differences be- 
tween the particle- and grid-based approach to solving Poisson's 
Equation, in an iV-body context. 

2 As we mention below, we actually analyse six simulations, but 
one of these is simply a higher-resolution version of one of the 
base simulations, and so is not entirely "independent" . We also 
include in our analysis a non cosmological simulation, see § 2.7. 
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Table 1. Summary of the properties for the simulation suite 



Simulation Name 


Code 




n 


ho 


n b 


thMM 


e 


Gas Resolution 


Feedback 






( M@) 








Gyr 


(kpc) 


( M @ ) 




S09 


RAMSES 


7.6xlO n 


0.3 


0.7 


0.045 


10.99 


0.4 


1.0 x 10 6 


kinetic 


G07(MW1) 


GASOLINE 


l.lxlO 12 


0.3 


0.7 


0.039 


10.89 


0.6 


8.0 x 10 5 


adiabatic 


B09(h277) 


GASOLINE 


7.1X10 11 


0.24 


0.77 


0.045 


10.41 


0.35 


1.6 x 10 4 


adiabatic 


B05(SGAL1) 


GCD+ 


5.0X10 11 


1 


0.5 


0.1 


9.74 


0.6 


1.0 x 10 6 


adiabatic 


A03 


GRAPESPH 


9.4xlO n 


0.3 


0.65 


0.045 


8.31 


0.5 


2.0 x 10 6 


kinetic 


H09 


RAMSES 


7.6X10 11 


0.3 


0.7 


0.045 


10.99 


0.2 


1.3 x 10 5 


adiabatic 


R08 


GASOLINE 


l.OxlO 12 


0.3 


0.7 


0.039 


N/A 


0.05 


0.2 x 10 5 


adiabatic 



into numerical simulations: one technique is to artificially 
turn off radiative cooling in the area where the SNell explo- 
sion occurs for a timescale long enough to allow the blast 
wave to expand. We call this type of feedback "adiabatic 
feedback". The second approach is to directly inject kinetic 
energy into the surrounding gas; we refer to this as "ki- 
netic feedback". Two of our simulations use kinetic feed- 
back, while the rest use an adiabatic approach. Pure "ther- 
mal feedback" is used in the case of Type la supernovae 
(SNela), where the longer lifetimes of the progenitors (rela- 
tive to SNell) means that the energy is not released into the 
same high density regions from which the stellar particles 
formed (and hence, the associated energy is not radiated 
away as efficiently as for the case of SNell). 

Each of the simulations employ star formation recipes 
which are similar; stars can form only from gas above a cer- 
tain density threshold. Since cosmological simulations typi- 
cally lack resolution below a few hundred parsecs, this sets 
a maximum density that the simulations can resolve on the 
order of 0.1 cm -3 ; all of the simulations here adopted this 
star formation threshold. We shall discuss the impact of this 
threshold selection later in the paper. 

The main properties of each of our simulations are pre- 
sented below, and summarised in Table 1. 



2.1 S09_YCosm_AMR_RAMSES 

We ran a high-resolution fully-cosmological (YCosm) disc 
simulation to redshift zero using the adap t ive m esh refine- 
ment (AMR)-based code RAMSES l|Tevssieri l|2002h l. The su- 
pernova feedback (SN) is modelled by directly injecting ki- 
netic energy into the surrounding gas - i.e., kinetic feedback. 

The simulation (S09, hereafter) was run within a "con- 
cordance" cosmology framework, with ilo = 0.3, ho = 
0.7, Qb = 0.045, and S1a = 0.7. Preliminary a nalysi s 
for this simulation was presented in iGibson et alJ (2009), 
while its optical proper t ies we re categorised extensively by 
ISanchez-Blazquez et al.l (2009). The simulated disc had its 
last major merger (LMM, defined as having total mass ratio 
of 1:3 or higher) at a redshift of z = 2.6, i.e. tLB = 10.99 Gyr 
(where tLB = lookback time), however, interactions with 
smaller satellites still occur at lower redshifts. We discuss the 
LMM later in the paper. Its final virial mass is 7.6x10" M 
at z — 0. 



2.2 G07_MWl_YCosm_SPH_GASOLINE 

T his galaxy is the s imula tion denoted as MW1 in the work 
of lGovernato et all (|2007l ). and is referred to as G07(MW1), 
hereafter. The code used for this fully-cosmological (YCosm) 
simulation is th e smoothed par t icle h ydrodynamics (SPH) 
code GASOLINE jWadslev et all lj2004T l). The SN feedback 
mechanism uses an adiabatic feedback approach where cool- 
ing was stopped artificially to allow blast waves from SNe to 
expand and heat the surrounding interstellar medium (ISM). 
In all the simulations run with GASOLINE presented here, 
40% of the SNe energy is coupled to the surrounding gas. 
Such a prescription results in a decrease in the amount of 
gas cooling early in the galaxy's formation, reducing the loss 
of angular momentum resulting from the merging of dense 
stellar systems. 

The simulation employed was run within a concordance 
cosmology with fio = 0.3, ho = 0.7, fib = 0.039, and fiA = 

0. 7; the last major merger was at redshift z=2.5, i.e. ti,B = 
10.89 Gyr, with several late minor interactions thereafter. 
The final virial mass is l.lxlO 12 M . 

2.3 B09_h277_YCosm_SPH_GASOLINE 

This simulation was also run with GASOLINE and was previ- 
ously studied in lBrooks et al.l (|2009h . and is referred to as 
B09(h277), hereafter. The simulation was run in a concor- 
dance cosmology, with fio = 0.24, ho = 0.77, fib = 0.045, 
and fi A = 0.76. The redshift of its LMM is also at z=2.5, 

1. e. t LB = 10.41 Gyr, but unlike the case for G07(MW1), 
this simulation experiences no mergers or accretion events 
since z«0.7 or tLB = 5.8 Gyr. The force resolution is also 
somewhat higher than for G07(MW1). 

2.4 B05_SGALl_NCosm_SPH_GCD+ 

This simulation was taken from the work of iBrook et al.l 
l|2005h . and is re ferred to as B05(SGAL1 ), hereafter. The 
SPH code GCD+ ijKawata fc Gibson! l|2003h ) was employed, 
although this particular run was not fully cosmological 
(NCosm). Semi-cosmological models, like B05(SGAL1), con- 
sist of an isolated sphere of dark matter and gas instead of 
a large cosmological volume. Small-scale fluctuations are su- 
perimposed on the sphere to allow for local collapse and sub- 
sequent star formation. Solid-body rotation is also applied 
to the sphere to incorporate the effects of longer wavelength 
fluctuations that a semi-cosmological model does not other- 
wise account for. Feedback from SNell was assumed to be 
adiabatic, with cooling turned off in the surrounding gas. 
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The cosmological framework in which B05(SGAL1) was 
run is quite different from the simulations discussed thus 
far; specifically, it used fio = 1, ho = 0.5, flf, = 0.1, and 
JIa = 0. While using the currently favoured ACDM frame- 
work would have a significant impact upon simulations of 
large-scale structure formation from Gaussian random noise 
initial conditions, it has been shown that, within the context 
of single galaxy formation models s uch as B05(SGAL1) , the 
resulting differences are negligible (|Brook et alJ (120051 1). In 
terms of its merger history, B05(SGAL1) is not dramatically 
different from B09(h277), in the sense of their being little or 
no merger activity since redshift z~0.5 or t^B = 5.9 Gyr. 

2.5 A03_YCosm_SPH_GRAPESPH 

This simul ation (hereaf t er refer red to as A03) was first pre- 
sented in lAbadi et al.l (|2003al ). It is a fully-cosmological 
(YCosm) Milky Way-like d i sc ga laxy, simulated with the 
GRAPESPH code l|Steinmetd (|l996h l. Feedback is predomi- 
nantly thermal, with 5% of supernova energy converted into 
kinetic feedback and injected into the surrounding gas par- 
ticles. A flat ACDM cosmology was assumed, with S7o = 0.3, 
ho = 0.65, fl;, = 0.045, and Qa = 0.7. Its final virial mass 
was 9.4XI0 11 Mq and its last major merger occurred at 
z = 1, i.e. tLB = 8.31 Gyr, although a number of minor 
interactions occur thereafter. 

2.6 H09_YCosm_AMR_RAMSES 

This simulation (H09, hereafter) traces the same halo as 
S09_YCosm_AMR_RAMSES described in § 2.1, but run with 
a higher spatial resolution, and employing a different feed- 
back mechanism for SNell. Instead of using kinetic feedback 
as in S09, it relies on an adiabatic feedback scheme. 

2.7 R08_NCosm_SPH_GASOLINE 

This s imulation (R08, hereafter) is taken from lRoskar et all 
(2008). It is an isolated Milky Way-type disc galaxy 
(l.OxlO 12 Mq) with solid body rotation added, similar to 
the B05_SGAll_NCosm_SPH_GCD+ simulation; where it 
differs from the latter is that R08 does not incorporate 
small scale density fluctuations. This means that this iso- 
lated simulation experiences no merger or accretion events. 
It is evolved for 10 Gyr using the same GASOLINE code as 
G07(MW1) and B09(h277), but at extremely high spatial 
resolution (softening length of 50 pc). 



3 RESULTS 

For observed solar neighbourhood stars, it is well established 
that there is a relationship between their velocity dispersion 
and age. We refer to the increase of the dispersion with time 
as disc heating, where the relationship for the solar neigh- 
bourhood indicates that the older disc stars are kinemati- 
cally "hotter" than its younger counterparts. Examining the 
dynamics of stars as a function of time therefore contains 
valuable information about the heating processes - driven by 
some combination of secular and satellite merger-driven phe- 
nomena. We focus on the vertical heating (a w , perpendicu- 
lar to the plane of the galaxy) as this out-of-plane heating is 



more susceptible to mergers/interactions. In-plane heating 
(ij u and o~ v ) is more sensitive to spiral wave and bar-driven 
heating which we do not consider in our studyQ 

In what follows, we first examine the velocity-dispersion 
age relation, similar to the manner by which observers study 
the same relation within the Milky Way, but now for stars 
within the simulated discs at 2=0. However, since the disper- 
sion of stars at z=0 does not provide direct information as 
to the velocity dispersion of the population at birth, we ex- 
tend our analysis to study the time evolution of the heating 
of sub-populations of disc stars. Specifically, we will attempt 
to ascertain whether stars which are kinematically hot to- 
day were born "hot" or were born "cold and heated" (by 
whatever means), thereafter. 

3.1 Age- Velocity Dispersion Relation 

We first examine the velocity dispersion perpendicular to the 
plane of the galaxy (a w ) for all stars at z = within the "lo- 
cal" disc, which we define as 7<R xy <9 kpc and \z\ < 1 kpc, 
as a function of age (Fig 1)0 With such a selection func- 
tion, all the simulations show little in the way of evidence 
for stellar heating for young stars with ages ~1 to ~3 Gyr, 
as well as much higher dispersions for old stars (consistent 
with observations of the Milky Way). 

We should point out that for the semi-cosmological sim- 
ulation, B05(SGAL1), we used a slightly larger radial cut of 
4<R xy <8 kpc and \z\ < 1 kpc due to the smaller number 
of star particles in this particular run. We found that in 
the smaller region of 7<R xy <9 kpc, there were not enough 
stars for some age bins to measure an accurate velocity dis- 
persion. We were therefore forced to use a larger region for 
this simulation. However, we also tested our results in all 
the simulations for three different radial cuts to ensure that 
the trends remain the same independent of radius. 

The S09 and H09 simulations show little sign of heating 
for stars younger than ~3 Gyr; for older stars, several dis- 
crete "jumps" in velocity dispersion can be discerned. Sim- 
ilar trends and small discrete jumps are also seen in the 
G07(MW1) simulation, in addition to a discontinuity in the 
dispersion for stars of age ~8 Gyr. The latter can be traced 
to a period of enhanced merger activity early in the galaxy's 
evolution, just prior to the establishment of its stable disc. 

Like all the simulations, A03 also shows dispersion 
trends consistent with the signature of continual/secular 
heating at later times. Stars older than ~9 Gyr, in partic- 
ular, possess a significantly large velocity dispersion. These 
high dispersions are a signature of the so-called "angular mo- 
mentum problem" mentioned in§ 1, which results in the for- 
mation of an overly-dominant spheroid compared to obser- 
vations. The spheroid component is dominated by old stars 
which are ultimately the responsible agents in the produc- 
tion of the high dispersions seen to the right-hand side of 
Fig 1. 

3 Further, Scabrokc & Gilmorc (2007) showed that dynamical 
streams can contaminate the local in-plane velocity distributions, 
which can complicate and compromise the comparison with simu- 
lated in-plane velocity distribution functions which do not capture 
adequately structure on that scale. 

4 In the sense that "young stars" have "small ages", in Fig 1, and 
"old stars" have "large ages". 
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Figure 1. Age- velocity dispersion relation in the vertical direc- 
tion for "local" disc stars in our suite of simulations. 



The absence of significant heating seen in all simulations 
for stellar ages of ~l-3 Gyr extends to somewhat older stars 
(up to ~6 Gyr in age), for both B09(h277) and B05(SGAL1). 
For older stars, discrete jumps in the dispersion, superim- 
posed upon a continual heating profile, are evident. The 
longer period during which stars show little heating is remi- 
niscent of the lQuillen fc Garnettl (|200lh interpretation of ex- 
tent observations that this reflects "saturation" in the thin 
disc's kinematic heating. The older, kinematically hotter, 
stars in these simulations have been suggested to be a sig- 
nature of the thick disc (e.g.. iBrook et al. (|2004h ). 

Broadly speaking, while there is a continuum of heat- 
ing "profiles" on display in Fig 1, at one end of the spec- 
trum, several of the simulations (e.g., S09, H09, G07(MW1)) 
show a temporal heating profile which becomes apparent at 
younger ages (~3 Gyr) relative to several at the opposite 
end of the spectrum (e.g., B05(SGAL1), B09(h277)) which 
only begin to show evidence of significant heating in older 
stars (>6 Gyr). If we associate these relatively flat periods 
at late times with the thin disc, then none of the simula- 
tions have thin discs as old as that of the Milky Way, which 
is considered to be between 8-10 Gyr old. 

One of the natural consequences of the merging which 
occurs within the hierarchical clustering paradigm is a de- 
gree of kinematical heating. As such, we set out to exam- 
ine the merging histories for each of the simulations, to see 
whether they shed light on the characteristic heating pro- 
files and discrete "jumps" seen in the age-velocity disper- 
sion plane (Fig 1); these merging histories will be discussed 
shortly. 

One thing which is readily apparent from our analy- 
sis is that all of the simulations show much higher velocity 
dispersions for the older stars in the disc, consistent with 
the behaviour seen in the Milky Way. That said, there is 
also a consistent offset, in the sense that all of the simulated 
discs are substantially hotter than that of the Milky Way. 
Part of this discrepancy relates to the fundamental prob- 
lem alluded to in § 1, specifically, that gas cools efficiently 
allowing baryons to collapse early during the merging pro- 
cess of galaxy formation, resulting in unrealistically large 
spheroidal components. These old spheroidal stars can have 



a significant impact on the derived dispersions, in the sense 
of "contaminating" what one would like to be a "pure" disc 
sample. In other words, rather than measuring the disper- 
sion of disc stars, which reflects the observational case, one 
is instead probing the additional impact that the dispersion 
of the spheroid stars have upon the sample. This is problem- 
atic, at some level, for all of the simulations - it is reflected 
in the very high dispersions seen at large ages in Fig 1. In 
order to make a fair comparison with disc stars from the 
Milky Way, we clean (in a very straightforward manner) our 
sample of these spheroidal contaminants by selecting in-situ 
stars, i.e. those that are born in the central galaxy. These in- 
situ stars are identified as those which form anywhere within 
the central galaxy, while stars that end in the central galaxy 
at z = but were born within a satellite or substructure are 
called accreted . We then derive the velocity dispersions of 
in-situ stars in our disc defined region, 7<R xy <9 kpc and 
\z\ < 1 kpc, at z = 0. 

By selecting in-situ stars we are examining the heating 
of disc stars. Whether this results in forming the thick disc 
or merely the old, hot thin disc, is left open to interpreta- 
tion. Of the thick disc formation mechanisms proposed (see 
§ 1), the direct accretion of satellites scenario is thus not 
explicitly addressed in this study. In this analysis we merely 
assume that a rotationally supported disc forms in-situ and 
can be born relatively hot or cold, and then may be heated 
by a number of processes. Further, recent simulation results 
have shown that some in - situ st ars will form part of the 
stellar halo l|Zolotov et al ] <|201Ch . and thus may affect our 
dispersion results. However, these stars are in the halo, with 
too few in our defined disc region to affect the dispersion-age 
plots presented here. 

We plot in Fig 2 the age-dispersion relation of these in- 
situ stars from the simulations along wi th three sets of ob- 
servat ional data for the Milky Way disc: Quillen fc Garnettl 
a combined set from ISoubiran fc Girardl (12005 



and 



Soubiran et al.l ll2008h , and Holmberg et al 



iQuillen fc Garnettl (l2001n" use a sample o f 189 nearby F- and 



G-dwarfs from Edvardsson et all (Il993l ); from their result- 
ing cr-age relation, they suggest that the Milky Way disc has 
been relatively quiescent with little heating for stars with 
ages between 3 and 9 Gyr, with stars older than that having 
been subject to an abrupt heatin g event. The second set o f 
obs ervational points is tak en from Soubiran fc Girardl (|2005T ) 
and lSoubiran et all (2008). We have merged these two cata- 
logues, in order to include a larger number of ol d disc stars in 
the sample. We note that th is dat a includes thelReddy et al.l 
(|2003h . lBeiisbv et all j|2003h and lBensbv et alj (2004) sam- 
ples, which target the thick disc specifically by using a kine- 
matic selection, and this may be the reason that their old 
stars are hotter than in the lHolmberg et al.l (|2008|) samples. 
Their analysis is consistent with that of Quillen fc Garnettl 
(2001), where the age-velocity relation of the thin disc is 
characterised by the saturation of the vertical dispersion at 
~25 km/s a t ages ^4—5 Gyr . The final set of observations 
is that from iHolmberg et al.l (|2008l ) ; they present a sample 
of F- and G-dwarfs from t he Geneva-Copenhagen Survey 
of the solar neighbourhood l|Nordstrom et al] (120041 ) . GCS) 
suggestive of a scenario consistent with continual heating of 
the local disc throughout its entire lifetime. 

By only considering in-situ disc stars, we have elim- 
inated a significant fraction of the high dispersion old 
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Figure 2. Vertical age- velocity dispersion relation for in-situ disc 
stars with in our suite of simulatio n s, compared with observa - 
tions fromlQuillen fc Garnetj J200lh: I Soubiran k. Girardl J2005f) ; 
ISoubiran et al. | teOOSft ; and lHolmberg et al. | ||2008| '). 



a force resolution of 50 pc, and it also exists in other 
high resolution isolated disc s i mulations in the li t eratur e 
(e.g.. iKazantzidis et all (120081 1: iKazantzidis et all (|2009h : 



IStewart et alj (|2009l )). In Fig 2, one can view our high- 
est resolution simulation (R08), as well as our lowest 
(B05(SGAL1)); if numerical heating was the main agent of 
the observed offset between the simulations and observa- 
tions, one might expect the lowest resolution simulation to 
be (kinematically) the hottest. This is not the case though 
and, in fact, B05(SGAL1) has the lowest resolution and is 
the coldest in the sample. 

Another important aspect to consider is the effect of 
secular heating; R08 has sufficient resolution to account 
for heating from internal processes such as from spiral 
arms. As the simulation is isolated and therefore removed 
from a cosmological context, the observed heating profile in 
this simulation must be secular due to spir al arms directly 
heatin g stars as well as causing migration ijLoebman et al.l 
(2010)). For the R08 simulation, these internal heating 
processes alone are e nough to match the observations of 
iHolmberg etU] (|200Sft . 



spheroidal components' contaminants; this is reflected in 
the ~20— 30% decrease in a w for stellar ages in excess at 
~7— 8 Gyr. We have only included four of the simulations 
in Fig 2, although the results that follow apply to the entire 
suite. The overall trend in Fig 2 matches that of Fig 1, in 
the sense that a range of heating "profiles" are seen, with 
both continual and discrete events being evident. 

We also compared the velocity dispersions in Fig 2 in 
three different regions for our highest resolution simulations 
(so H09, R08 and B09(h277), with selected regions being 
4 < R xy < 8 kpc, 6 < R xy < 9 kpc, and 7 < R xy < 9 
kpc). These simulations have enough stars to examine much 
smaller volume cuts. We found that the trends in velocity 
dispersion remain the same, independent of radial cut, with 
quantitative differences with Fig 2 being insignificant. 

It might be argued that the three simulations 
with the (relatively) oldest disc component (B09(h277), 
B05(SGAL1), R08) - also, those with relatively flat cr-age 
relations for ages up to ~6 Gyr - are a somewhat better 
reflection of the relation inferred from the observations. We 
will show in § 3.3 that the merger history of these systems 
is a primary driver in the establishment of this relationship, 
and examine the role played by numerical effects. 

An interesting observation from Fig 2 is that there is 
a significant offset, even when including just in-situ stars, 
with all simulations compared with observations at any time 
in the galaxy's evolution, in the sense that the stars in 
the simulations are hotter than the observed stars at all 
times (with the exception of the semi-cosmological simula- 
tion, B05(SGAL1)). This offset is particularly high when 
looking at old stars but is also significant for young stars. 
Several possibilities might be responsible for driving such a 
discrepancy: (i) numerical heating due to limited force reso- 
lution, (ii) the treatment of heating and cooling within the 
ISM of the simulations, and (iii) the adopted low star for- 
mation threshold. 

The issue of numerical heating will be addressed at 
length in § 4; here, we si mply note that the offset also 
exists in the simulation of [Roskar et al.l (|2009h , which has 



3.2 Are Stars Born Hot or Heated Subsequently? 

In this section, we aim to answer several questions that 
emerged from the above kinematical analysis of z = 
stars: were the kinematically hot, old, stars in Figs 1 and 
2 born with these high velocity dispersions, or were they 
born "cold" and heated subsequently? If the latter, then 
what is the source of this heating? To answer these, we ex- 
amine the kinematics of disc stars at the time of their birth 
for different epochs of a galaxy's formation. We do this by 
selecting disc stars born in the "disc" region, A<R xy <& kpc 
and z<l kpc, at the time of their birth, using a fairly ar- 
bitrary time "slice" of 200 Myrs - i.e., we are deriving the 
velocity dispersion of young disc stars in each simulation at 
various epochs. We tested our results with different radial 
cuts and age range, and found that our results and conclu- 
sions are not sensitive to the used values. The slightly larger 
radial slice used in this section allows us to obtain a larger 
sample of stars to derive their dispersions. 

Figure 3 shows the derived velocity dispersions for 
young stars at different times throughout the respective sim- 
ulations' evolution. Each of the orthogonal components of a 
are highlighted, although as noted earlier, our analysis will 
concentrate solely upon a m . For clarification, stars born at 
early times are situated to the left of each panel in Fig 3, 
while stars born more recently are located towards the right 
- i.e., the abscissa now reflects "cosmic time" rather than 
"stellar age" (as was employed in Figs 1 and 2). 

For the S09 simulation (top left panel of Fig 3) , all disc 
stars, independent of time, are born cold with low vertical 
velocity dispersions of a w ~ 30 kms -1 , on average. There is 
a slight increase in the dispersion for stars with formation 
times between ~7 and 10 Gyr, where the dispersion increases 
by ~25%. This epoch corresponds to a period of enhanced 
minor merger activity, during which the ISM is heated kine- 
matically relative to the adjoining quiescent phases. 

For the G07(MW1) simulation (top middle panel of 
Fig 3), stars are born on average with vertical velocity dis- 
persions between a w = 20 and a w — 30 kms^ 1 , except for 
the period between ~3 and 5 Gyr. During this time there are 
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Figure 3. Radial, tangential, and vertical velocity dispersion components (<r u , cr v , and cr,,,, respectively) for young stars (<200 Myrs) 
at the time of their birth for different epochs, within the defined disc annulus 4 < R xy < 8 kpc and \z\ < 1 kpc. 



several minor mergers with satellites which result in these 
stars being born with velocity dispersions roughly twice that 
of the adjoining phases (a m = 60 kms -1 ). It is also impor- 
tant to note that these mergers produce a short-lived warp 
at z = 2 - i.e., at t = 3.2 Gyr. The stars that we detected 
in the disc during this period were located within this warp 
region. Because of their potential to dominate over in-plane 
stars at only a few scalelengths, stars in the warp should be 
treated carefully, particularly in the case of studying their 
kinematics, as they ca n result in an appa rent increase in 
the velocity dispersion (|Roskar et al.l l|201(t )). These "warp" 
stars are kinematically "disturbed" and born with higher 
u w . This is a very similar trend to that seen in the H09 
simulation (bottom right panel of Fig 3) where stars, on 
average, tend to have dispersions between a w = 20 and 
<T W = 30 kms" 1 at the time of birth, but there is a pe- 
riod between ~3 and 5 Gyr, again, where this dispersion 
doubles to about a w = 60 kms" 1 . As for G07(MW1), this 
period coincides with an epoch of enhanced satellite interac- 
tion with the main galaxy, although in G07(MW1) the warp 
is the primary cause of the high velocity dispersion during 
this period and not the minor mergers. 

A distinct trend is noticed for the B09(h277), 
B05(SGAL1), and A03 simulations. Stars born at late times 
(over the past ~6 Gyr) are born cold, with velocity disper- 
sions between 10 and 20 kms -1 , while stars born prior to this 
are born hot (with vertical velocity dispersion of ~70 km/s). 
It is tempting to interpret this as the signature of separate 
thin and thick discs, where the thick disc is composed of 
older star s which were born h otter than the younger, colder, 
thin disc ijBrook et all {200l)). 

The high velocity dispersion measured at early times 
in A03, i.e. from ~2 to 6 Gyr, is due to the numer- 
ous merge r events that this si mulations undergoes at early 
times fsee lAbadi et al.l (|2003al 'll. The feedback mechanism 
is not particularly effective and therefore the satellites that 
merge with the main galaxy contain a large stellar compo- 



nent which affect the high velocity dispersions derived. The 
merger activity is largely over by ~6 Gyr and the disc is al- 
lowed to settle and form. One can interpret the low velocity 
dispersions determined in Fig 3 for stars from t = 6 Gyr as 
signature of the formation of such disc. 

Having identified the velocity dispersions of stars at 
birth, we now wish to determine whether they maintain the 
self-same dispersion as they age - i.e., are these stars being 
heated with time? We do this by selecting the same "young" 
stars at a particular time and then tracing them forward in 
time, in order to quantify the degree of evolution in the ve- 
locity dispersion of these ensembles of stars. This is shown 
in Fig 4, where the subsequent heating of the stars at each 
epoch is represented by the coloured curves. Because we are 
interested in stars born in the disc of the galaxy we nec- 
essarily choose epochs after the disc has formed, with the 
exception of A03 where the disc forms much later compared 
to the other simulations. For each galaxy this time can vary 
depending upon the time of the last major merger. We there- 
fore do not look at stars beyond z ~ 2.5 because, in general, 
the discs in these galaxies have not yet formed. 

Looking first at the S09 simulation, the stars born at 
t = 2.5 Gyr have an initial velocity dispersion of a w = 30 
kms -1 , increasing to a w — 70 kms -1 over the subsequent 
~2 Gyr. This behaviour is qualitatively repeated for all stars 
born (and tracked) in the first ~6 Gyr: i.e., stars are born 
relatively cold but rapidly heat to more than double their 
initial velocity dispersion within ~1 Gyr, before the heating 
begins to "saturate" , while stars born over the past ~6 Gyr 
(while also born "cold") heat much more gradually. Indeed, 
stars born over the last ~4 Gyr experience essentially no 
kinematic heating. 

G07(MW1) presents a qualitatively similar heating pro- 
file to that of S09, in the sense of (a) older stars expe- 
riencing a doubling of their vertical velocity dispersion in 
the first few Gyr after birth, before the heating saturates, 
and (b) younger stars experiencing little, if any, kinemati- 
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Figure 4. Velocity dispersion perpendicular to the plane (a w )ol young stars (ages <200 Myrs) at various epochs (represented by different 
colours, and different starting times), traced forward in time to quantify their temporal heating profile. 



cal heating. One subtle difference between G07(MW1) and 
S09, though, is that older stars in the former are born rela- 
tively hot compared with their younger counterparts, while 
in S09, all stars, independent of birth epoch, are born with 
essentially the same vertical velocity dispersion. 

A03 is similar to G07(MW1) in the sense that older 
stars are born relatively hot compared with the younger 
stars in the disc. The difference in this simulation compared 
to both S09 and G07(MW1) is that the younger stars also 
experience heating, where their dispersions double in the 
first few Gyr. It can be seen that even stars born at t = 13 
Gyr are being heated with time (see red curve in bottom 
middle panel of Fig 4) . 

Both B09(h277) and B05(SGAL1) present somewhat 
different heating profiles. With the exception of the first 
epoch in both (coinciding with the time of the last major 
merger in both cases), the vertical velocity dispersion is es- 
sentially invariant - i.e., older stars (which are born hotter 
than their contemporary counterparts, as in G07(MW1)) 
and younger stars maintain their birth velocity dispersion 
for the lifetime of the simulation. 

The next question which needs addressing is this: what 
are the responsible agents driving the assorted heating pro- 
files seen in Fig 4? What is the quantitative relationship to 
their respective merger histories? Is the effect of warps play- 
ing an important role? Are numerical effects plaguing the 
analysis? We address these over the following sub-sections. 

3.3 The Effects of Mergers 

As noted above, our efforts have concentrated upon the out- 
of-plane heating within the simulations, due to its stronger 
sensitivity to mergers and interactions. It is crucial to derive 
and quantify the merger histories of our simulations, in order 
to link the observed heating trends with the interactions 
they have experienced during their evolution. 

We have already seen indirect signatures of mergers in 



the above analysis. Stars born with hotter kinematics during 
early stages of a galaxy's evolution (Fig 3), as well as the 
dramatic heating of stars born at early times over fairly short 
timescales (Fig 4), can be related to the last major merger 
(LMM) of each galaxy (see Table 1 for the lookback time, 
thMM, at which the LMM occurred for each simulation). 
Major mergers have a total mass ratio of 3:1 or higher. 

In what follows, we examine minor mergers of satellites 
with mass 4% the mass of the disc at the time of the merger, 
back to z ~ 2. Such merg ers are able to disturb disc struc- 
ture (|Quinn et al.l (jl993)). We are, unfortunately, limited 
by time resolution due to the available number of outputs 
for each simulation. We are thus not able to trace directly 
the trajectories of the satellites and determine whether they 
penetrate the disc, or the number of close passages which 
occur prior to the final coalescence. We restrict our analy- 
sis to satellites which have contributed stars to the inner 10 
kpc of the central galaxy by z — 0, indicating that these 
satellites have interacted with the disc. 

The last major merger in the S09 simulation occurs at 
lookback time of thMM = 10.99 Gyr, which corresponds to 
a time t = 2.02 Gyr in Fig 3. It has a mass ratio of 3:1 
(Af„ir=l-lxl0 u and Af ao t=1.7xlQ 10 ). This major merger 
heats the disc stars significantly as can be seen in the black 
line in the top left panel in Fig 3. After the LMM there 
are several minor baryonic mergers, the most noteworthy of 
which occurs between redshifts z=0.8 and z=0.7 (a time cor- 
responding to ~6.3 Gyr in Fig 3). This minor merger has a 
mass ratio of 8:1 {M vir =5.7x 10 11 and M sat =6xl0 10 ). Stars 
born during this period are somewhat hotter kinematically- 
speaking, relative to those born before and after (Fig 3). In 
addition, stars born in the preceding ~3 Gyr to this merger 
appear to have been subject to rapid heating (see the yel- 
low and cyan lines in upper left panel of Fig 4) . Additional 
(less significant in terms of mass) mergers occur at redshift 
z = 1.75 — 1.44 - i.e., t — 4 — 4.5 with a mass ratio of 16:1. 
The effects of these mergers are not obvious in our plots, 
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although the impact of the former likely plays a role in the 
heating seen between times 3 and 4 Gyr in the upper left 
panel of Fig 4. This simulation undergoes minor mergers up 
to redshift z — 0. 

G07(MW1) undergoes its LMM at lookback time 
thMM = 10.89 Gyr, which corresponds to time t = 2.5 Gyr 
in Fig 3. This major merger has a mass ratio of 3:1 
(M„ ir =9xl0 10 and M sat =2.5x 10 10 ), where the heating 
caused by this merger can be seen in the increase in the 
dispersion for stars during this period (see black line in top 
middle panel from Fig 3). It undergoes several minor merg- 
ers at time t = 3.2 Gyr with mass ratios of 10:1 and 15:1, 
the effects of which can also been seen in Fig 3 and Fig 4. 
Hereafter, it undergoes minor interactions the last one oc- 
curring at t = 7.68 Gyr. These interactions are apparent 
near t = 8 Gyr in Figs 3 and 4, although their heating "im- 
pact" is not particularly obvious in Fig 3 - i.e., the stars 
born at t=8 Gyr (upper middle panel of Fig 3) are not kine- 
matically hotter than those born within ~2 Gyr of these 
mergers; similarly, the impact on the heating of these stars 
is not particularly dramatic (upper middle panel of Fig 4). 
The effect that the LMM (at i~2.5 Gyr) has in heating 
the ISM of G07(MW1) is felt over the subsequent ~3 Gyr 
(Fig 3), and impacts upon the temporal heating profiles of 
stars born during this period as well (Fig 4). G07(MW1) 
hosts fairly significant warps during and after these peri- 
ods of merger-driven activity, the impact of which has been 
noted previously (§ 3.2). 

For B09(h277) there is a clear distinction between stars 
born at early times and those formed at later times, which 
can be ascribed to the simulation's merger history. The 
LMM in this simulation occurs at a formation time of 
thMM = 10.41 Gyr, corresponding to time t — 3.2 Gyr 
in Fig 3 and has a mass ratio of 3:1 (M„i r =6x 10 10 and 
M sat =2x 10 10 ). During the period between t=2 and 3.4 Gyr, 
a significant number of both major and minor interactions 
take place, with a final baryonic interaction at t ~ 3.4 Gyr 
(with a mass ratio of 100:1). This period of merger activity 
maps directly onto the time during which the ISM is signifi- 
cantly hotter (upper right panel of Fig 3). The complete lack 
of major or minor baryonic mergers subsequent to this point 
is reflected in the absence of detectable temporal heating in 
stars born since £~4 Gyr. 

We have already discussed the similarity between the 
heating profiles of B05(SGAL1) and B09(h277), where there 
is a clear distinction between stars born at early epochs and 
those born later. B05(SGAL1) has its LMM at a forma- 
tion time of thMM = 9.74 Gyr or t~3.3 Gyr in Fig 3 and 
has a mass ratio of 3:1 (M vir =6xl0 10 and M sat =2x 10 10 .) 
It experiences only one minor 10:1 (Mi« r =6xl0 10 and 
M sat =2xl0 9 ) baryonic interaction at t~6 Gyr after the 
LMM, the signature of which is not readily seen in Figs 3 or 
4. 

A03 undergoes the last major merger with mass ratio 
3:1 at t~6 Gyr and it lasts for ~1 Gyr. The ISM is hotter 
during this LMM phase, as evidenced in the higher a w at 
time t~6.5 Gyr in the lower middle panel of Fig 3. At times 
earlier than Gyr, there are many major merger events 
as can be seen by the large velocity dispersions measured for 
these stars. It undergoes it last merger event at z~0.74, i.e. 
t=7.5 Gyr in Fig 3 and 4 and has completely merged with 
the disc of the main galaxy by z~0.48, i.e t~9.2 Gyr, with 



a mass ratio of 45:1 (M vir =9x 10 11 and M 3at =2x 10 10 ). See 
lAbadi et all (|2003bh for details of this satellite. The heat- 
ing profile of stars formed more recently in A03 (i.e., those 
formed within the final 3—4 Gyr of the simulation) differs 
from those of the other simulations, in sense that even these 
recently-formed stars within A03 experience significant heat- 
ing. This can be associated to a companion satellite that 
survives at z=0. It first appears within a radius of 15 kpc at 
2=0.33, i.e., £=10.52 Gyr, with a mass of M sat =5.8x 10 9 . 

We separate our simulations into two groups, those 
which have interactions at low redshift after the thin disc has 
formed and those that show no major or minor interactions 
since redshift z~l. S09 (and H09), G07(MW1), and A03, 
undergo later minor merger interactions and therefore ex- 
hibit more evidence of continual, later, heating. Conversely, 
B09(h277), B05(SGAL1), and R08 (not shown) experience 
no later minor merger activity once the disc has formed, 
and therefore do not show jumps in the heating over short 
timescales in their respective discs at later times, which is 
associated with such merger events in the other simulations. 
If we combine this information with what we deduce from 
looking at the age-velocity dispersion plane in Fig 1 and 
Fig 2, we can conclude that in order to obtain a thin disc 
consistent with observations, the simulated galaxy must ex- 
perience no interactions at late times (at least, since 2~1). 

3.4 The Central Concentration of the Satellites 

The effect of heating that accreted satellites have on the disc 
is dependent on the mass distribution of the satellite, in the 
sense that the accretion of more massive, and more concen- 
trated satellites, wi ll cause a higher degree of heating (e.g. 
IVelazauez fc White! ll 19991 ) ). Simulations produce rotation 
curves that rise rapidly in t he inner regions wi th a central 
peak before dropping off (e.g. lMaver et al.l ([2008)). However, 
observations of dwarf galaxies have shown that their rota- 
tion curves rise linearly in the central regions. Presumably, 
accreted satellites should have mass distributions which are 
similar to local galaxies. The more concentrated satellites in 
the simulations are related to the "angular momentum prob- 
lem" , where the baryons are deficient in angular momentum 
and produce overly concentrated stellar bulges. This chal- 
lenge for cold dark matter cosmology is beyond the scope of 
this paper, but we note that several mechanisms have been 
proposed to reso lve the discrepancy b etween theory and 
observ a tion (e.g. Navarro et all (I 19961): Mashchenko et all 
l|2008h : IScannapieco et al l ^OOg )). iGovernato et al l (|20ich 
showed that resolution which is high enough to form lo- 
cal star formation within an inhomogeneous ISM, will drive 
large scale supernova outflows and decrease the central mass 
concentration, producing simulated dwarfs which have a 
mass distribution that matches observed galaxies. The reso- 
lution required to create such dwarfs is not achieved in any 
simulation of a Milky Way mass galaxy in our study, or in- 
deed in the literature. 

We looked at the rotation curves of satellites in three 
of the simulations discussed in this paper, chosen at redshift 
z ~ 2, each with a dynamical mass of ~10 10 M and they 
showed peaked rotation curves indicative of an excess of cen- 
tral material. We have shown that the major source of disc 
heating in our suite of simulations is due to the interaction 
and accretion of satellite galaxies with the disc. The high 
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central mass concentration of our satellites may be causing 
these effects to be exaggerated compared to the effect of 
real satellites, particularly if such satellites do indeed have 
"cored" r ather than cuspy central density profiles and no 
bulge (e.g. lOh et al.l l|2010t) ). This effect is perhaps the most 
important caveat to our work; future, increased, resolution 
which results in more r ealistic dwarfs (akin to those seen in 
iGovernato et al.l (|2010l 11 may reduce the heating rates seen 
in the current suite of cosmological simulations. 



4 THE EFFECTS OF RESOLUTION AND 
STAR FORMATION RECIPES 

It is important to determine whether numerical heating is 
influencing our results. Two-body heating can have a dra- 
matic effect on the increase in the k inetic e n ergy of a kine- 
matically cold rotating stellar disc ifM avcr (2004])) and is, 
therefore, an important factor to take into consideration in 
our stu dy. Such numerical e ffects are dependent upon reso- 
lution l|Moore et all (|l996h : ISteinmetz fc White] l|l997h ). 

In our sample of simulations we have a variety of res- 
olutions. S09, G07(MW1), B05(SGA11), and A03 have rel- 
atively low spatial resolution - between 400 and 600 pc - 
while B09(h277) and H09 have somewhat higher resolution 
(~300 pc). The isolated disc from R08 has a much higher 
resolution (50 pc). If the heating we see was dominated by 
numerical effects, one might expect a particularly large effect 
in the lowest resolution simulation: B05(SGAL1). In fact, 
this is the coolest of all the simulations studied. Further, we 
have shown that B09(h277), B05(SGAL1), and R08 present 
similar trends in the heating of their disc stars during their 
quiescent period of evolution at low redshift, showing little 
stellar heating, despite having vastly different resolutions. 
Of the simulations which show significant recent merger ac- 
tivity, H09 is the highest resolution, yet it shows heating at 
low redshift in agreement with the lower resolution simula- 
tions which have similar merging histories (S09, G07(MW1), 
and A03). These trends appear to indicate that numerical 
heating is not the main driver of the inferred heating profiles. 

However, the importance of the spectre of numerical 
heating means that one must proceed cautiously and exam- 
ine the issue more quantitatively. In cosmological simula- 
tions, resolution dependence is more complicated than the 
case where isolated discs are used as the initial conditions. 
There is numerical heating due to gravitational softening, 
but on the other hand when we go to higher resolution we 
resolve more substructure, creating more heating. Another 
problem is that low resolution substructures tend to be (ar - 
tificially) more concentrat ed (|van den Bosch et all l^OOlT) ; 
iBarnes fc Hernauistl (|l996h '). meaning that the heating ef- 
fects of their interactions may be exaggerated. We discussed 
the central concentrations of the satellites in these cosmo- 
logical simulations in § 3.4. 

In order to explore possible numerical heating effects, 
we therefore examined a set of isolated disc galaxies. The 
initial conditions were created as in iKazantzidis et aD 
(2008), to which the reader is referred for details, and 
were run using GASOLINE. The isolated galaxies that 
we re-simulate compris es an exponentia l stellar disc, a 
Hernquist model bulg e dHernauist |l990|)'), and an NFW 
dark matter profile l|Navarro et alj (|l997h ). The total 



mass of the galaxy is 10 12 M©, similar to that of the 
Milky Way, with a disc gas fraction of 10%. To form a 
rotationally-supported disc we impart angular momentum 
to the gas component corresponding to a spin parameter 
of A=0.04. We evolve each simulation for 1 Gyr, after 
allowing 0.2 Gyr for the system to relax. The disc has 
an initial Toomre stability parameter equal to Q = 2.2, 
which means it is stable against any local nonaxisymmetric 
instabilities. The main differences between the isolated runs 
are summarised in Table 2. We run three simulations at 
different resolutions - high (ISO _HR_LT_GASOLINE), 
medium (ISO^IRXT.GASOLINE), and low 
(ISO_LR_LT_GASOLINE) - and we employ the same 
star formation threshold (0.1 cm -3 ) used in the cosmo- 
logical simulations analysed in this study. The purpose of 
this is to see the effects resolution might have on heating 
stars in simulations. We then run another high resolution 
simulation, but employ a much higher star formation 
threshold (100 cm" 3 : ISO_HR_HT_GASOLINE), more 
comparable to the densities associated with star formation, 
observationally. 

In Fig 5, we plot the age-dispersion relation for stars 
within our set of isolated Milky Way disc galaxies. While 
there is a not surprising resolution-dependency in the verti- 
cal velocity dispersions of the stars at birth (ranging from 
~16 km/s, to ~21 km/s, to ~29 km/s, for the HR, MR, 
and LR runs, respectively), there is little, if any, evidence 
for any significant heating within any of the simulations, ir- 
respective of their different resolutions. We will return to 
the differences that resolution has upon the stellar velocity 
dispersions at birth, at the end of this section. 

The differences in the velocity dispersions between the 
isolated runs seen in Fig 5 could be due to the effects of 
feedback in these simulations. In order to determine how 
important this effect may be when determining the disper- 
sion in Fig 5, we compare the star formation rates (SFR) 
in all the isolated runs. We find that the dispersion is not 
greatly affected by the feedback in the sense that, for ex- 
ample, changes in the SFR (by factors of between 2 and 4) 
did not change the dispersions. Further, comparison between 
the low and high threshold dispersion at times when they 
had equivalent SFRs, we find the same offsets as indicated 
in Fig 5. 

The highest resolution, isolated simulation, like the case 
of R08, is particularly interesting in the context of this sec- 
tion. Both these have significantly higher resolution than the 
cosmological simulations. ISO JHR_LT .GASOLINE and R08 
are both isolated from a cosmological context, and so no 
heating from satellites occurs, yet these simulations have a 
vertical velocity dispersion "floor" of ~ 15— 20 kms" 1 , similar 
to the "floor" in dispersions seen in the fully cosmological 
simulations during their respective quiescent periods. R08 
uses the same star formation and supernova feedback physics 
as G07(MW1) and B09(h277); this dispersion "floor" is tied 
directly to the implementation of ISM physics within the 
code. Such physics is difficult to capture in cosmological 
simulations, as it is multi-scale, going from kpc-scale, where 
most of the gas is ionised, to pc-scale, where most of the gas 
is molecular. Cosmological simulations, like the ones anal- 
ysed here, follow the formation of galaxies in a volume of 
at least several tens of Mpc, because aspects of structure 
formation require that the large scale gravitational field is 
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Table 2. Summary of the Properties for the Isolated Simulations 



Simulation Name 


Star Mass 


Gas Mass 


Dark Mass 






Star Formation Threshold 




(M ) 


(M ) 


(M Q ) 


(kpc) 


(kpc) 


(cm -3 ) 


ISO_HR_HT_GASOLINE 


1.73x10 s 


2.56xl0 4 


1.35x10 s 


0.1 


0.17 


100 


ISO_HR_LT .GASOLINE 


1.73X10 3 


2.56xl0 4 


1.35x10 s 


0.1 


0.17 


0.1 


ISO_MR_LT .GASOLINE 


1.38X10 4 


2.05xl0 5 


1.08x10 s 


0.2 


0.34 


0.1 


ISO.LRXT.GASOLINE 


1.04x10 s 


1.63xl0 6 


8.68x10 s 


0.4 


0.68 


0.1 



a Dark matter softening length 
^ Gas softening length 



properly modelled. Related to this, our inability to resolve 
locally collapsing high density regions means that we aver- 
age star formation over large columns, using a low density 
threshold for star formation (0.1 cm -3 in the simulations 
analysed here). Yet star formation is observed to occur in 
regions where gas has cooled to regions of significantly higher 
density. The low star formation threshold means that, within 
the simulations, gas may be forming stars in regions which 
remain relativ ely hot. 

Recently, iGovernato et al] (|2010l ) showed that with a 
spatial resolution of about ~100 pc, gas could be allowed 
to collapse to densities more representative of the average 
density observed in star forming giant molecular clouds. Us- 
ing a star formation density threshold of 100 cm -3 , they 
successfully simulated the f irst bulgeles s disc galaxy (see 
IGovernato et alj (|2010T ) and iBrook et all (|2010h ). However, 
if the density threshold is too high for the resolution of the 
simulation, the galaxy comes out too compact. We imple- 
ment these High Threshold (HT) recipes within a high res- 
olution isolated simulation (ISO JHR_HT_GASOLINE) , and 
overplot the age- velocity dispersion relationship (cyan trian- 
gles) in Fig 5. Two striking features are immediately appar- 
ent: (i) the dispersion is much lower than for all the simula- 
tions which used a low star formation density threshold, even 
when using the same high resolution. The difference is far 
greater than the difference which was caused by resolution^ 
(ii) very little heating occurs, even at these very low disper- 
sions. This shows very clearly that numerical heating is not 
affecting our cosmological simulations. Rather, a dispersion 
floor is created by the inability of gas to cool sufficiently be- 
fore forming stars when the star formation density threshold 
is set at a lower level (0.1 cm -3 ). 

The resolution-dependent differences in the dispersions 
of the isolated simulations are due to the differences in the 
degree to which the gas is able to cool before forming stars. 
This is supported by three pieces of information: (i) no heat- 
ing is apparent at any of the three significantly different res- 
olutions, with the age-dispersion relation remaining flat; (ii) 
stars form with lower dispersion at higher resolution; (iii) 
even at the very low dispersion levels of the High Thresh- 
old simulation (~5 kms -1 ), heating was negligible. These 
conclusions are supported by the fact that the temperature 
of the gas from which stars are born increases, as resolu- 



5 sec Pilkington et al. 2010, in prep, for a detailed analysis of the 
ISM velocity dispersion as a function of star formation threshold 
and resolution. 
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Figure 5. Age-velocity dispersion relation in the vertical direc- 
tion for disc stars (4 < R xy < 8 kpc and \z\ < 1 kpc) for four 
isolated Milky Way-scale simulations. 



tion decreases, with the average temperature being 7300 K, 
6500 K, 5900 K and 400 K, respectively, for the LR, MR, 
HR, and HRJHT isolated galaxies. 



5 SUMMARY 

We have analysed the kinematics of disc stars in a suite of 
Milky Way-scale simulations which were run with different 
hydrodynamical cosmological codes and at different resolu- 
tions. Some were run using the Smooth Particle Hydrody- 
namics (SPH) approach whereas others used the Adaptive 
Mesh Refinement (AMR) method. This is the first paper to 
compare cosmological disc galaxies run with these two dif- 
ferent techniques. No differences in the analysed kinematic 
properties of the simulated galaxies were found to be de- 
pendent on the approach used for the implementation of 
gas hydrodynamics. 

First, we analysed the velocity dispersion of all disc 
stars as a function of age at z = 0, comparing with anal- 
ogous observations of the Milky Way's disc. An overall off- 
set exists, in the sense that all the simulated galaxies are 
hotter than the Milky Way's disc. This was shown to be 
driven in part by resolution and star formation threshold ef- 
fects, although the latter is much more efficient at reducing 
the dispersion "floor" (Fig 5). We provide evidence that the 
dominant contributor is the low density threshold for star 
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formation, which has been routinely implemented in sim- 
ulations of Milky Way-scale galaxies, although we should 
point out that we only show this for non-cosmological sim- 
ulations. This low density threshold means that stars are 
formed from unphysically high temperature gas, creating a 
dispersion "floor". Indeed, our lowest resolution simulation 
has the least amount of kinematical heating. Future cos- 
mological simulations with sufficient resolution to resolve 
the mean density of g iant molecular clouds (akin to the 
iGovernato et all (|2010T ) dwarf galaxy simulations) will be 
a critical step forward in this work. 

Despite this dispersion "floor" in our simulations, some 
interesting heating trends are found. Two of our simulations 
(B09(h277) and B05(SG AL1)) are in better agreem ent with 
interpretations made bv lQuillen fc Gam ctt (200l]), where a 
saturated disc is present for young stars up to t = 6 Gyr and 
discrete jumps seen in the dispersion for older disc stars. The 
other simulations (S09, G07(MW1), A03, and R08) seem 
to be in better agreement with the disc undergoing con- 
tinuou s heating, consistent with analysis of iHolmberg et aD 
(2008), although the rate of heating in the simulations re- 
mains higher than that observed in nature. 

We then proceeded to study the heating of these stars 
as a function of time; starting from the point at which the 
final disc was stable (z~2), we derived the dispersion of 
stars at the time of their birth and how those coeval en- 
sembles evolved with time. We found that whereas in some 
simulations stars are born cold in the disc and are heated 
(S09, G07(MW1), A03, and H09), either numerically or 
due to a physical process, in other simulations (B09(h277), 
B05(SGAL1) and R08) the stars maintain essentially the 
same dispersion as they possessed at birth. Further, in some 
simulations, stars are born with high dispersio ns - i.e., they 
are bo rn "hot" . This ca n be due to in t eracti ons ijBrook et al.l 
I^OOJ ) 1 ) and/or warps (|Roskar et alJ (|2010h ). Turbulence in 
the interstellar medium (ISM) not related to mergers could 
also be a cause of stars being born with large velocity disper- 
sion. Recent observations of high redshift discs indicate that 
internal proc esses are a po s sible ca use of the observ e d tur - 
bulent ISM (|Genzel et alJ l^OOSh l. iBournaud et~aH l|2009h 
compared simulations formed internally in unstable gas-rich, 
clumpy discs with simulations of merger induced disc thick- 
ening, and found that thick discs formed internally are a bet- 
ter match to observed high redshift discs. Mechanisms such 
as cold flows and supernova feedback are currently being dis- 
cussed - in addition to mergers- as possi ble causes of t he tur - 
bulent ISM in high red s hift systems ( Bju^grtetaL l|2010T l; 
iForster Schreiber et all l^OlcMCeverino et al.l (|2010l )V 

Within the favoured cosmological paradigm of hierar- 
chical clustering, merging and accretion of satellites onto 
host galaxies is fundamental. Our goal has been to examine 
the effects that these mergers might have upon the heat- 
ing of disc stars. We find a clear relationship when looking 
at the heating profiles between those simulations that have 
late mergers and those that heat significantly. Four simula- 
tions (S09, G07(MW1), A03 and H09) have minor mergers 
at low redshift, and we map these interactions onto the in- 
creases seen in the velocity dispersion of their disc stars. 
The other three simulations (B09(h277), B05(SGAL1), and 
R08), which have no interactions over the past ~7 Gyr, show 
little heating in the disc with time. We note that R08 is an 
isolated simulation which has no satellites, and hence no in- 



teractions, yet has sufficient heating due to spiral arms and 
migration to match observed heating rates of the Milky Way. 
The suite of cosmological simulations do not have the ability 
to resolve these secular effects, nor heating due to molecular 
clouds. In these simulations, it is only in the quiescent period 
since the last accretion events that heating is low enough to 
match the Milky Way's thin disc. None has a thin disc older 
than ~ 6 Gyr, indicating that it would be difficult to gain 
a thin disc as old as some estimates for the Milky Way thin 
disc within the current cold dark matter paradigm. A caveat 
of our study is the overly concentrated mass distributions of 
our satellites, meaning that resolution of this persistent "old 
thin disc" problem may come from improved modelling of 
baryonic physics coupled with increased resolution. 
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